 rm(list = ls())


 library(car)  # needed for recode function
 library(Zelig)
 library(Amelia)

#setwd("")

load("d_RUS_FIRMS_2010_MI.RData")


library(sandwich) 
library(lmtest) 


######  LAWYERS2



b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(lawyers2~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results

coef3.lawyers2<-combined.results$q.mi[1]
se3.lawyers2<-combined.results$se.mi[1]

######  COURTS2


b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(courts2~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.courts2<-combined.results$q.mi[1]
se3.courts2<-combined.results$se.mi[1]


######  LAWENF2


b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(lawenf2~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.lawenf2<-combined.results$q.mi[1]
se3.lawenf2<-combined.results$se.mi[1]


######  BUR2


b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(bur2~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.bur2<-combined.results$q.mi[1]
se3.bur2<-combined.results$se.mi[1]


######  COURTS_INF2


b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(courts_inf2~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.courts_inf2<-combined.results$q.mi[1]
se3.courts_inf2<-combined.results$se.mi[1]


######  LAWENF_INF2




b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(lawenf_inf2~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.lawenf_inf2<-combined.results$q.mi[1]
se3.lawenf_inf2<-combined.results$se.mi[1]


######  BUR_INF2



b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(bur_inf2~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.bur_inf2<-combined.results$q.mi[1]
se3.bur_inf2<-combined.results$se.mi[1]

######  PSS2



b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(pss2~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.pss2<-combined.results$q.mi[1]
se3.pss2<-combined.results$se.mi[1]



######  PSA2


b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(psa2~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.psa2<-combined.results$q.mi[1]
se3.psa2<-combined.results$se.mi[1]



######  MAFIA2




b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(mafia2~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.mafia2<-combined.results$q.mi[1]
se3.mafia2<-combined.results$se.mi[1]



######  LAWYERS



b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(lawyers~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.lawyers<-combined.results$q.mi[1]
se3.lawyers<-combined.results$se.mi[1]


######  COURTS


b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(courts~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.courts<-combined.results$q.mi[1]
se3.courts<-combined.results$se.mi[1]


######  LAWENF



b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(lawenf~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.lawenf<-combined.results$q.mi[1]
se3.lawenf<-combined.results$se.mi[1]



######  BUR



b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(bur~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.bur<-combined.results$q.mi[1]
se3.bur<-combined.results$se.mi[1]


######  COURTS_INF


b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(courts_inf~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.courts_inf<-combined.results$q.mi[1]
se3.courts_inf<-combined.results$se.mi[1]


######  LAWENF_INF


b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(lawenf_inf~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.lawenf_inf<-combined.results$q.mi[1]
se3.lawenf_inf<-combined.results$se.mi[1]



######  BUR_INF




b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(bur_inf~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.bur_inf<-combined.results$q.mi[1]
se3.bur_inf<-combined.results$se.mi[1]

######  PSS



b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(pss~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.pss<-combined.results$q.mi[1]
se3.pss<-combined.results$se.mi[1]



######  PSA



b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(psa~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.psa<-combined.results$q.mi[1]
se3.psa<-combined.results$se.mi[1]


######  MAFIA


b.out<-NULL
se.out<-NULL

for(i in 1:RF$m) 
{
JGM.out <- lm(mafia~privatized*consolidated+tax1_dich_r+other_firms2_r+factor(size_dummies)+firm_age+foreign_own+gov_own+factor(growth)+BA+age+gender+leg_ed+factor(job_descript)+rts_viol+litigated+city2+factor(sector3), data = RF$imputations[[i]])
VC<-vcovHC(JGM.out, type="HC1")
intercoeff<-JGM.out$coeff[3]+JGM.out$coeff[41]
intersd<-sqrt(VC[3,3]+VC[41,41]+2*VC[3,41])
b.out <- rbind(b.out, intercoeff)
se.out <- rbind(se.out, intersd)
}
combined.results <- mi.meld(q = b.out, se = se.out)
combined.results


coef3.mafia<-combined.results$q.mi[1]
se3.mafia<-combined.results$se.mi[1]

###############
#FIGURE 3 CODE
###############
# will produce PDF output

#pdf("fig_3", height =6, width = 9) #open pdf device


layout(matrix(c(3,1,2), 1,3),widths=c(1,3,2.4))  ## allows for three graphs, in matrix function is "plot1" through "plot 3"; then commands for 1 rows, 3 columns; c(3,1,2) puts 3rd plot first; widths sets relative width of columns


##PLOT 1

#Create Vectors for coefs and standard errors for each model, and variable names
  

coef.vec.1 <- c(coef3.mafia2,coef3.psa2,coef3.pss2,coef3.courts_inf2,coef3.lawenf_inf2,coef3.bur_inf2,coef3.lawyers2,coef3.courts2  ,coef3.lawenf2,coef3.bur2)
 se.vec.1 <- c(se3.mafia2,se3.psa2,se3.pss2,se3.courts_inf2,se3.lawenf_inf2,se3.bur_inf2,se3.lawyers2,se3.courts2,se3.lawenf2,se3.bur2)

 
var.names <- c("criminal racket", "private security \nagency ","internal security \nservice ", "courts \n(with connections)","law enforcement \n(unofficial)","bureaucrats \n(unofficial)","lawyers \n(out of court)","courts", "law enforcement \n(official)", "bureaucrats \n(official)")

y.axis <- length(var.names):1#create indicator for y.axis, descending so that R orders vars from top to bottom on y-axis
adjust <- .18 #create object that we will use to adjust points and lines up and down to distinguish between models
    
par(mar=c(3,8,7,1), lheight = .8)#set margins for regression plot, c(bottom, left, top, right)
plot(coef.vec.1, y.axis, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = 1.3, #plot model 1 coefs using black points (pch = 19, default = black), adding the "adjust amount" to the y.axis indicator to move points up

xlim = c(-2.1, 2.1), xaxs = "r", #set limits of x-axis so that they include mins and maxs of 
ylim = c(min(y.axis)-.2, max(y.axis)+.2),yaxs="r", main = "",cex.main=1)
   title("Property Dispute",line=5.8,cex.main=1.7)
axis(1,at = seq(-4,4, by = 1), label = seq(-4,4, by = 1), mgp = c(0,1,0), cex.axis = 1.2)#add x-axis and labels; c(title,labels,line) "pretty" creates a sequence of  equally spaced nice values that cover the range of the values in 'x'-- in this case, integers
axis(2, at = y.axis, label = var.names, las = 1, tick = T, cex.axis =1.2)#add y-axis and labels; las = 1 makes labels perpendicular to y-axis
#axis(3,pretty(coef.vec.1, 3))#same as x-axis, but on top axis 
abline(h = y.axis, lty = 2, lwd = .5, col = "light grey")#draw light dotted line at each variable for dotplot effect
box(bty="l")#draw box around plot
segments(coef.vec.1-qnorm(.975)*se.vec.1, y.axis, coef.vec.1+qnorm(.975)*se.vec.1, y.axis, lwd =  1.3)#draw lines connecting 95% confidence intervals
#for tick marks indicating 90% ci's use following 2 lines:
#segments(coef.vec.1-qnorm(.95)*se.vec.1, y.axis -.035, coef.vec.1-qnorm(.95)*se.vec.1, y.axis +.035, lwd = 1.1)#draw vertical tick marks at 90% 
#segments(coef.vec.1+qnorm(.95)*se.vec.1, y.axis -.035, coef.vec.1+qnorm(.95)*se.vec.1, y.axis+.035, lwd = 1.1)   #confidence intervals
abline(v=0, lty = 2) # draw dotted line through 0 for reference line for null significance hypothesis testing
points(coef.vec.1, y.axis, pch = 21, cex = 1.3, bg = "white" ) #add point estimates for 2nd model; pch = 21 uses for overlay points, and "white" for white color


arrows(.2,10.8,1.5,10.8, length=.07, lty=1, lwd=.7, xpd=TRUE)
arrows(-.2,10.8,-1.5,10.8, length=.07, lty=1, lwd=.7, xpd=TRUE)

text(.2,11.45, "consolidated owners \nmore likely to use than \n unconsolidated owners",xpd=TRUE,adj=0,cex=1.2)
text(-.2,11.45, "consolidated owners \nless likely to use than \n unconsolidated owners",xpd=TRUE,adj=1,cex=1.2)

                
#### PLOT 2


coef.vec.1 <- c(coef3.mafia,coef3.psa,coef3.pss,coef3.courts_inf,coef3.lawenf_inf,coef3.bur_inf,coef3.lawyers,coef3.courts  ,coef3.lawenf,coef3.bur)
 se.vec.1 <- c(se3.mafia,se3.psa,se3.pss,se3.courts_inf,se3.lawenf_inf,se3.bur_inf,se3.lawyers,se3.courts,se3.lawenf,se3.bur)

y.axis <- length(var.names):1#create indicator for y.axis, descending so that R orders vars from top to bottom on y-axis
adjust <- .18 #create object that we will use to adjust points and lines up and down to distinguish between models
    
par(mar=c(3,2,7,2), lheight = .8)#set margins for regression plot, c(bottom, left, top, right)
plot(coef.vec.1, y.axis, type = "p", axes = F, xlab = "", ylab = "", pch = 19, cex = 1.3, #plot model 1 coefs using black points (pch = 19, default = black), adding the "adjust amount" to the y.axis indicator to move points up
 
 xlim = c(-2.1, 2.1), xaxs = "r", #set limits of x-axis so that they include mins and maxs of 
 ylim = c(min(y.axis)-.2, max(y.axis)+.2),yaxs="r", main = "",cex.main=1)
   title("Debt Dispute",line=5.8,cex.main=1.7)
axis(1,at = seq(-4,4, by = 1), label = seq(-4,4, by = 1), mgp = c(0,1,0), cex.axis = 1.2)#add x-axis and labels; c(title,labels,line) "pretty" creates a sequence of  equally spaced nice values that cover the range of the values in 'x'-- in this case, integers
#axis(2, at = y.axis, label = "", las = 1, tick = T, cex.axis =1)#add y-axis and labels; las = 1 makes labels perpendicular to y-axis
#axis(3,pretty(coef.vec.1, 3))#same as x-axis, but on top axis 
abline(h = y.axis, lty = 2, lwd = .5, col = "light grey")#draw light dotted line at each variable for dotplot effect
#box(bty="l")#draw box around plot
segments(coef.vec.1-qnorm(.975)*se.vec.1, y.axis, coef.vec.1+qnorm(.975)*se.vec.1, y.axis, lwd =  1.3)#draw lines connecting 95% confidence intervals
#for tick marks indicating 90% ci's use following 2 lines:
#segments(coef.vec.1-qnorm(.95)*se.vec.1, y.axis -.035, coef.vec.1-qnorm(.95)*se.vec.1, y.axis +.035, lwd = 1.1)#draw vertical tick marks at 90% 
#segments(coef.vec.1+qnorm(.95)*se.vec.1, y.axis -.035, coef.vec.1+qnorm(.95)*se.vec.1, y.axis +.035, lwd = 1.1)   #confidence intervals
abline(v=0, lty = 2) # draw dotted line through 0 for reference line for null significance hypothesis testing

points(coef.vec.1, y.axis, pch = 21, cex = 1.3, bg = "white" ) #add point estimates for 2nd model; pch = 21 uses for overlay points, and "white" for white color

arrows(.2,10.8,1.5,10.8, length=.07, lty=1, lwd=.7, xpd=TRUE)
arrows(-.2,10.8,-1.5,10.8, length=.07, lty=1, lwd=.7, xpd=TRUE)

text(.2,11.45, "consolidated owners \nmore likely to use than \n unconsolidated owners",xpd=TRUE,adj=0,cex=1.2)
text(-.2,11.45, "consolidated owners \nless likely to use than \n unconsolidated owners",xpd=TRUE,adj=1,cex=1.2)


##########
library(pBrackets)

par(mar=c(3,.5,7,1.5))#set margins for regression plot, c(bottom, left, top, right)

#left.side <- .45#use this to manipulate how far segments are from y-axis
    #note:  getting braces and text in proper place requires much trial and error

plot(seq(0,1,length=length(var.names)), y.axis, type = "n", axes = F,  xlab = "", ylab = "")#call empty plot using type="n"


brackets(1, 1.2, 1, 4, h = NULL, ticks = 0.5, curvature = 0.5, type = 1, col = 1, lwd = 1, lty = 1, xpd = TRUE)
text(.55, 1.6, labels=expression(paste("Legal Strategies")), adj=c(0,0),xpd=TRUE,srt=90, cex=1.2)


brackets(1, 5, 1, 7, h = NULL, ticks = 0.5, curvature = 0.5, type = 1, col = 1, lwd = 1, lty = 1, xpd = TRUE)
text(.55, 5.05, labels=expression(paste("Illegal Strategies")), adj=c(0,0),xpd=TRUE,srt=90,cex=1.2)
text(.72, 5.3, labels=expression(paste("(corruption)")), adj=c(0,0),xpd=TRUE,srt=90,cex=1.2)

brackets(1, 7.9, 1, 9.8, h = NULL, ticks = 0.5, curvature = 0.5, type = 1, col = 1, lwd = 1, lty = 1, xpd = TRUE)
text(.55, 7.9, labels=expression(paste("Illegal Strategies")), adj=c(0,0),xpd=TRUE,srt=90,cex=1.2)
text(.72, 8.3, labels=expression(paste("(violence)")), adj=c(0,0),xpd=TRUE,srt=90,cex=1.2)


#dev.off()





